## load packages

library(readxl)
library(ggplot2)
install.packages("margins")
library(margins)
library(dplyr)
library(dynlm)
install.packages("zoo")
library(zoo)
install.packages("ARDL")
library(ARDL)
install.packages("plm")
library(plm)
## MODEL 1
data <- read_excel("data_ardl.xlsx")
View(data)
names(data)[names(data) == "Public_it-1"] <- "public_previous"

ts_data <- ts(data$Laws_it, start = c(2002), end = c(2023), frequency = 1)  # Assuming 'Laws_it' is the 4th column
model_1_ardl <- lm(Laws_it ~ public_previous + Misery_index + Loc_elect + Gen_elect + Presidentialism + Electoral_Index, data = data)
summary(model_1_ardl)

model_poisson <- glm(Laws_it ~ Electoral_Index, family = poisson(link = "log"), data = data)
new_data <- data.frame(Electoral_Index = seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 100))
predicted_counts <- predict(model_poisson, newdata = new_data, type = "response")
predicted_counts_mean <- predict(model_poisson, newdata = new_data, type = "link", se.fit = TRUE)
confidence_interval <- exp(cbind(predicted_counts_mean$fit - 1.96 * predicted_counts_mean$se.fit, 
                                 predicted_counts_mean$fit + 1.96 * predicted_counts_mean$se.fit))
plot(new_data$Electoral_Index, predicted_counts, type = "l", col = "blue", lwd = 2,
     xlab = "Electoral Democracy Index", ylab = "Predicted Counts of Laws Enacted")




selected_points <- plot_data[seq(1, nrow(plot_data), length.out = 10), ]

ggplot(data = plot_data, aes(x = Electoral_Index, y = Predicted_Probability)) +
  geom_line(color = "black", linetype = "dashed") +  # Main trend line
  geom_ribbon(aes(ymin = Predicted_Probability - 1.96 * Standard_Error, 
                  ymax = Predicted_Probability + 1.96 * Standard_Error),  # Shadow bands
              fill = "grey", alpha = 0.3) +  # Specify shadow color and transparency
  geom_point(data = selected_points, aes(x = Electoral_Index, y = Predicted_Probability), size = 2) +  # Selected points
  geom_errorbar(data = selected_points, aes(x = Electoral_Index, 
                                            ymin = Predicted_Probability - 1.96 * Standard_Error, 
                                            ymax = Predicted_Probability + 1.96 * Standard_Error),  # Vertical lines for confidence intervals
                width = 0.05, color = "black") +  # Specify line width and color
  xlab("Electoral Democracy Index") +
  ylab("Predicted Probability of Laws Enacted") +
  ggtitle("Predictive Marginal Effect (Model 1)")

### MODEL 2

library(dynlm)
ts_data <- ts(data$Laws_it, start = c(2002), end = c(2023), frequency = 1)  # Assuming 'Laws_it' is the 4th column

model_2_ardl <- dynlm(Laws_it ~ L(public_previous, 1) + Misery_index + Loc_elect + Gen_elect + Electoral_Index, data = data)
summary(model_2_ardl)
## model 3 

ts_data <- ts(data$Laws_it, start = c(2002), end = c(2023), frequency = 1)
View(data)
library(plm)
install.packages("marginaleffects")
library(marginaleffects)
# Assuming 'entity' represents entities and 'time' represents time periods

## model 3
random_effects_model3 <- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + democ, data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model3)


## model 2
random_effects_model2 <- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + Electoral_Index, data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model2)

## model 1

random_effects_model1<- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + Presidentialism +Electoral_Index , data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model1)


## plot
# Fit the Poisson regression model
model_poisson <- glm(Laws_it ~ Electoral_Index, family = poisson(link = "log"), data = data)

# Generate new data with sequence of values for Electoral_Index
new_data <- data.frame(Electoral_Index = seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 100))

# Predict counts of Laws_it based on new data
predicted_counts <- predict(model_poisson, newdata = new_data, type = "response")

# Predicted counts and confidence interval
predicted_counts_mean <- predict(model_poisson, newdata = new_data, type = "link", se.fit = TRUE)
confidence_interval <- exp(cbind(predicted_counts_mean$fit - 1.96 * predicted_counts_mean$se.fit, 
                                 predicted_counts_mean$fit + 1.96 * predicted_counts_mean$se.fit))

# Create a dataframe for plotting
plot_data <- data.frame(Electoral_Index = new_data$Electoral_Index,
                        Predicted_Probability = predicted_counts,
                        Standard_Error = predicted_counts_mean$se.fit)

# Extract 10 points for plotting
selected_points <- plot_data[seq(1, nrow(plot_data), length.out = 10), ]

# Plot the predictive marginal effect with selected points and vertical lines
ggplot(data = plot_data, aes(x = Electoral_Index, y = Predicted_Probability)) +
  geom_line(color = "black", linetype = "dashed") +  # Main trend line
  geom_ribbon(aes(ymin = Predicted_Probability - 1.96 * Standard_Error, 
                  ymax = Predicted_Probability + 1.96 * Standard_Error),  # Shadow bands
              fill = "grey", alpha = 0.3) +  # Specify shadow color and transparency
  geom_point(data = selected_points, aes(x = Electoral_Index, y = Predicted_Probability), size = 2) +  # Selected points
  geom_errorbar(data = selected_points, aes(x = Electoral_Index, 
                                            ymin = Predicted_Probability - 1.96 * Standard_Error, 
                                            ymax = Predicted_Probability + 1.96 * Standard_Error),  # Vertical lines for confidence intervals
                width = 0.05, color = "black") +  # Specify line width and color
  xlab("Electoral Democracy Index") +
  ylab("Predicted Probability of Laws Enacted") +
  ggtitle("Predictive Marginal Effect (Model 1)")




# Fit the random effects model 2
random_effects_model2 <- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + Electoral_Index, 
                             data = data, index = c("CAP_domain", "Years"), model = "random")

# Extract the coefficients from the model summary
coefficients <- coef(random_effects_model2)

# Extract the standard errors of the coefficients
se <- sqrt(diag(vcovHC(random_effects_model2)))

# Calculate the predicted probabilities and confidence intervals
predicted_probs <- plogis(coefficients["Electoral_Index"] + seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 10) * coefficients["Electoral_Index"])

# Create a dataframe for plotting
plot_data <- data.frame(Electoral_Index = seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 10),
                        Predicted_Probability = predicted_probs)

# Calculate the length of the extended vertical lines
vertical_line_length <- 0.9 * (max(predicted_probs) - min(predicted_probs))

# Plot the predictive marginal effect with 10 points and extended vertical lines
ggplot(data = plot_data, aes(x = Electoral_Index, y = Predicted_Probability)) +
  geom_line(color = "black", linetype = "dashed") +  # Main trend line
  geom_point(size = 2) +  # Add points
  geom_errorbar(aes(ymin = Predicted_Probability - vertical_line_length * se["Electoral_Index"], ymax = Predicted_Probability + vertical_line_length * se["Electoral_Index"]),  # Extended vertical lines
                width = 0.05, color = "black") +  # Specify line width and color
  xlab("Electoral Democracy Index") +
  ylab("Predicted Probability of Laws Enacted") +
  ggtitle("Predictive Marginal Effect (Model 2)")


##### third model plot

# Fit the random effects model
random_effects_model3 <- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + democ, 
                             data = data, index = c("CAP_domain", "Years"), model = "random")

# Extract the coefficients from the model summary
coefficients <- coef(random_effects_model3)

# Extract the standard errors of the coefficients
se <- sqrt(diag(vcovHC(random_effects_model3)))

# Calculate the predicted probabilities and confidence intervals
predicted_probs <- plogis(coefficients["democ"] + seq(min(data$democ), max(data$democ), length.out = 10) * coefficients["democ"])

# Create a dataframe for plotting
plot_data <- data.frame(democ = seq(min(data$democ), max(data$democ), length.out = 10),
                        Predicted_Probability = predicted_probs)

# Calculate the length of the extended vertical lines
vertical_line_length <- 0.4 * (max(predicted_probs) - min(predicted_probs))

# Plot the predictive marginal effect with 10 points and extended vertical lines
ggplot(data = plot_data, aes(x = democ, y = Predicted_Probability)) +
  geom_line(color = "black", linetype = "dashed") +  # Main trend line
  geom_point(size = 2) +  # Add points
  geom_errorbar(aes(ymin = Predicted_Probability - vertical_line_length * se["democ"], ymax = Predicted_Probability + vertical_line_length * se["democ"]),  # Extended vertical lines
                width = 0.05, color = "black") +  # Specify line width and color
  xlab("Interaction Term (Presidentailism*Electoral Index)") +
  ylab("Predicted Probability of Laws Enacted") +
  ggtitle("Predictive Marginal Effect (Model 3)")





#### model 1 again

# Fit the random effects model
random_effects_model1 <- plm(Laws_it_panel ~ public_previous + Misery_index + Loc_elect + Gen_elect + Presidentialism + Electoral_Index, 
                             data = data, index = c("CAP_domain", "Years"), model = "random")

# Extract the coefficients from the model summary
coefficients <- coef(random_effects_model1)

# Extract the standard errors of the coefficients
se <- sqrt(diag(vcovHC(random_effects_model1)))

# Calculate the predicted probabilities and confidence intervals
predicted_probs <- plogis(coefficients["Electoral_Index"] + seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 10) * coefficients["Electoral_Index"])

# Create a dataframe for plotting
plot_data <- data.frame(Electoral_Index = seq(min(data$Electoral_Index), max(data$Electoral_Index), length.out = 10),
                        Predicted_Probability = predicted_probs)

# Calculate the length of the extended vertical lines
vertical_line_length <- 0.4 * (max(predicted_probs) - min(predicted_probs))

# Plot the predictive marginal effect with 10 points and extended vertical lines
ggplot(data = plot_data, aes(x = Electoral_Index, y = Predicted_Probability)) +
  geom_line(color = "black", linetype = "dashed") +  # Main trend line
  geom_point(size = 2) +  # Add points
  geom_errorbar(aes(ymin = Predicted_Probability - vertical_line_length * se["Electoral_Index"], ymax = Predicted_Probability + vertical_line_length * se["Electoral_Index"]),  # Extended vertical lines
                width = 0.05, color = "black") +  # Specify line width and color
  xlab("Electoral Democracy Index") +
  ylab("Predicted Probability of Laws Enacted") +
  ggtitle("Predictive Marginal Effect (Model 1)")




#### some other results if you use lagged value of public opinion $$$
library(stats)
library(dplyr)

random_effects_model1 <- plm(Laws_it_panel ~ lag(public_previous, 1) + Misery_index + Loc_elect + Gen_elect + Presidentialism + Electoral_Index, 
                             data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model1)


random_effects_model2 <- plm(Laws_it_panel ~ lag(public_previous, 1) + Misery_index + Loc_elect + Gen_elect + Electoral_Index, data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model2)


random_effects_model3 <- plm(Laws_it_panel ~ lag(public_previous, 1) + Misery_index + Loc_elect + Gen_elect + democ, data = data, index = c("CAP_domain", "Years"), model = "random")
summary(random_effects_model)



